#include "uhm.hxx"

#include "uhm/wrapper/fort/util.hxx"
#include "uhm/wrapper/fort/common.hxx"
#include "uhm/wrapper/fort/linal.hxx"

using namespace uhm;

#ifndef __ID_FORT__
#define __ID_FORT__ 1
#endif 

#define N __ID_FORT__

typedef ID_Node_H_<N>          Id_Nd;
typedef ID_Element_<1>         Id_El;

// Solvers

// -- LU 
F_void UHM_C2F(uhm_direct_lu_nopiv_initialize)(F_int *type, F_int *rhs, F_int *blk, 
                                               F_ptr *solver) { 
  Direct_Solver_Base_ *s 
    = new (std::nothrow) Direct_Solver_<Id_Nd,Id_El,LU_nopiv_Linal_>
    (*type, Mat_Base_::FULL, Mat_Base_::UNSYM, *rhs, *blk);
  UHM_ERROR(s != NULL, UHM_ERROR_SOLVER_NEW);
  solver->first  = (F_long)s;
  solver->second = UHM_FORT_LU_NOPIV;
}

F_void UHM_C2F(uhm_direct_lu_piv_initialize)(F_int *type, F_int *rhs, F_int *blk, 
                                             F_ptr *solver) { 
  Direct_Solver_Base_ *s
    = new (std::nothrow) Direct_Solver_<Id_Nd,Id_El,LU_piv_Linal_>
    (*type, Mat_Base_::FULL, Mat_Base_::UNSYM, *rhs, *blk);
  UHM_ERROR(s != NULL, UHM_ERROR_SOLVER_NEW);
  solver->first  = (F_long)s;
  solver->second = UHM_FORT_LU_PIV;
}

F_void UHM_C2F(uhm_direct_lu_incpiv_initialize)(F_int *type, F_int *rhs, F_int *blk, 
                                                F_ptr *solver) { 
  Direct_Solver_Base_ *s
    = new (std::nothrow) Direct_Solver_<Id_Nd,Id_El,LU_incpiv_Linal_>
    (*type, Mat_Base_::FULL, Mat_Base_::UNSYM, *rhs, *blk);
  UHM_ERROR(s != NULL, UHM_ERROR_SOLVER_NEW);
  solver->first  = (F_long)s;
  solver->second = UHM_FORT_LU_INCPIV;
}

// -- LDL
F_void UHM_C2F(uhm_direct_ldl_nopiv_initialize)(F_int *type, F_int *rhs, F_int *blk,
                                                F_ptr *solver) {
  Direct_Solver_Base_ *s
    = new (std::nothrow) Direct_Solver_<Id_Nd,Id_El,LDL_nopiv_Linal_>
    (*type, Mat_Base_::LOWER_TRIANGULAR, Mat_Base_::SYM, *rhs, *blk);
  UHM_ERROR(s != NULL, UHM_ERROR_SOLVER_NEW);
  solver->first  = (F_long)s;
  solver->second = UHM_FORT_LDL_NOPIV;
}

F_void UHM_C2F(uhm_direct_ldl_piv_initialize)(F_int *type, F_int *rhs, F_int *blk,
                                              F_ptr *solver) {
  Direct_Solver_Base_ *s
    = new (std::nothrow) Direct_Solver_<Id_Nd,Id_El,LDL_piv_Linal_>
    (*type, Mat_Base_::LOWER_TRIANGULAR, Mat_Base_::SYM, *rhs, *blk);
  UHM_ERROR(s != NULL, UHM_ERROR_SOLVER_NEW);
  solver->first  = (F_long)s;
  solver->second = UHM_FORT_LDL_PIV;
}

F_void UHM_C2F(uhm_direct_ldl_incpiv_initialize)(F_int *type, F_int *rhs, F_int *blk,
                                                 F_ptr *solver) {
  Direct_Solver_Base_ *s
    = new (std::nothrow) Direct_Solver_<Id_Nd,Id_El,LDL_incpiv_Linal_>
    (*type, Mat_Base_::LOWER_TRIANGULAR, Mat_Base_::SYM, *rhs, *blk);
  UHM_ERROR(s != NULL, UHM_ERROR_SOLVER_NEW);
  solver->first  = (F_long)s;
  solver->second = UHM_FORT_LDL_INCPIV;
}

// -- CHOL
F_void UHM_C2F(uhm_direct_chol_nopiv_initialize)(F_int *type, F_int *rhs, F_int *blk,
                                                 F_ptr *solver) {
  Direct_Solver_Base_ *s
    = new (std::nothrow) Direct_Solver_<Id_Nd,Id_El,CHOL_nopiv_Linal_>
    (*type, Mat_Base_::LOWER_TRIANGULAR, Mat_Base_::SYM, *rhs, *blk);
  UHM_ERROR(s != NULL, UHM_ERROR_SOLVER_NEW);
  solver->first  = (F_long)s;
  solver->second = UHM_FORT_CHOL_NOPIV;
}

F_void UHM_C2F(uhm_direct_chol_piv_initialize)(F_int *type, F_int *rhs, F_int *blk,
                                               F_ptr *solver) {
  Direct_Solver_Base_ *s
    = new (std::nothrow) Direct_Solver_<Id_Nd,Id_El,CHOL_piv_Linal_>
    (*type, Mat_Base_::LOWER_TRIANGULAR, Mat_Base_::SYM, *rhs, *blk);
  UHM_ERROR(s != NULL, UHM_ERROR_SOLVER_NEW);
  solver->first  = (F_long)s;
  solver->second = UHM_FORT_CHOL_PIV;
}

F_void UHM_C2F(uhm_direct_chol_incpiv_initialize)(F_int *type, F_int *rhs, F_int *blk,
                                                  F_ptr *solver) {
  Direct_Solver_Base_ *s
    = new (std::nothrow) Direct_Solver_<Id_Nd,Id_El,CHOL_incpiv_Linal_>
    (*type, Mat_Base_::LOWER_TRIANGULAR, Mat_Base_::SYM, *rhs, *blk);
  UHM_ERROR(s != NULL, UHM_ERROR_SOLVER_NEW);
  solver->first  = (F_long)s;
  solver->second = UHM_FORT_CHOL_INCPIV;
}

